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Aims. We present a newly developed code that allows simulations of optical observations of galaxy fields with a variety of instruments. The 
code incorporates gravitational Iensing eff'ects and is targetted at simulating Iensing by galaxy clusters. Our goal is to create the tools required 
for comparing theoretical expectations with observations to obtain a better understanding of how observational noise aff'ects Iensing applications 
such as mass estimates, studies on the internal properties of galaxy clusters and arc statistics. 

Methods. Starting from a set of input parameters, characterizing both the instruments and the observational conditions, the simulator provides 
a virtual observation of a patch of the sky. It includes several sources of noise such as photon-noise, sky background, seeing, and instrumental 
noise. Ray-tracing through simulated mass distributions accounts for gravitational Iensing. Source morphologies are realistically simulated 
based on shapelet decompositions of galaxy images retrieved from the GOODS-ACS archive. According to their morphological class, spectral- 
^ energy-distributions are assigned to the source galaxies in order to reproduce observations of each galaxy in arbitrary photometric bands. 

Results. We illustrate our techniques showing virtual observations of a galaxy-cluster core as it would be observed with the space telescope 
rS DUNE, which was recently proposed to ESA within its "Cosmic vision" programme. We analyze the simulated images using methods applicable 
^ I to real observations and measure the properties of gravitational arcs. In particular, we focus on the determination of their length, width and 
" " ' curvature radius. 

Conclusions. We find that arc properties strongly depend on several properties of the sources. In particular, our results show that compact, faint 
or low surface brightness galaxies that are barely detectable are more easily distorted as arcs with large length-to-width ratios. We conclude 
that realistic Iensing simulations can be obtained with the method proposed here. They will be essential for evaluating and improving analysis 
techniques currently used for cosmological interpretations of cluster Iensing. 

Key words, gravitational Iensing - galaxies: clusters, dark matter 



1. Introduction 

Gravitational Iensing of distant galaxies is one of the most pow- 
erful tools for investigating the content of the universe and for 
understanding how the matter is distributed within the cosmic 
structures. 

Being the most massive objects in the universe, galaxy clus- 
ters are the most powerful existing gravitational lenses. This 
makes them suitable for a large number of purposes. 

Clusters are able to magnify distant galaxies. Thus, they 
are natural gravitational telescopes, that can be used to investi- 



gate the first sources of light in the universe (lEllis et alJ l2QQll: 
Schaerer & PeiiQll20Q3l: iKneib et alJEooil: iLi et al.ll2QQ7 



Strong Iensing features, like gravitational arcs and multi- 
ple images are observed in the central, most dense re gions of 
the lenses dFort et al.lll988[ iLvnds & Petrosianlll989h . Hence, 
they can be used t o constrain the inner struc ture of clusters 



fe.g ILi et alJ I20Q7I: iC omerford et al. 2006; Q guri et alJ I2QQ3 



Sand et al.ll2Q04t iMeneghetti et al.ll2007LI2QQ5al) 



Gravitational arcs are also used in a variety of cosmologi- 
cal applications. Their frequency on the sky reflects the abun- 
dance, the concentration and the redshift distribution of mas- 
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sive lenses in the universe. Moreover, it has been shown that the 
pro bability of forming arcs is enhanced during cluster merg- 



ers 



dTorri et al.1 12004 iFedeli et al.1 12OQ6). Thus, arc statistics 



are a powerful tool for tracing structure formation and for con- 
straining the cosmological parameters dBartelmann et al.ll 19981: 
Li et al .l l2005l:lMeneghetti et allEoOSbl) . 



Moreover, since they form along critical lines, whose size 
scales with the source and lens redshifts, lensed images aris- 
ing from multiple sources at different redshift s can be used 
to constrain the geometry of the universe (Gols e et aD E002: 
[Meneghetti et al.ll2005a|) . and thus also cosmological parame- 
ters. 

In the outer regions of clusters, the distortions induced 
by the cluster potenti als are weaker. Nevertheless, u sing suit- 
able algorithms (e.g. iBartelmann & Schneidej|200ll and ref- 
erences therein), the distortion field can be inverted to obtain 
two dimensional s urface density maps of the the lenses (see 



Clowe et al.ll2004l: IStarck et al.l 120061: IPaulin-Henriksson et al.. 
2OO7I for some recent examples). Strong and weak tensing can 



also be used jointly to improve the mass reconstrution in a 
non-parametric vv ^ay (|Bradac et al.ll2004l: ICacciato et al.ll2006l: 
Diego et al.ll2007l) . 



Although the potential power of gravitational tensing for 
the above mentioned applications has been already demon- 
strated, diff'erent implementations of methods adopted for ex- 
tracting the desired information from the tensing data often 
give contradictory results. An example is given by the so called 
"arc statistics problem", i.e. the claimed inconsistency between 
the observed and the predicted number of highly distorted grav- 
itational arcs expec ted to be seen on the sky in the standard 



ACDM cosmology Je^telmann et al]ll998l; IPalalet a l. 2005; 
Wambsganss et al.l l2004 iLuppino eTaP Il999l: [( 



IGladders et al.l 



20031: IZaritskv & GonzalezI 120031: iLi et al.1 hOOSh. Moreover, 
testing a particular analysis technique and/or comparing the re- 
sults to the theoretical expectations is difficult because obser- 
vational noise int roduces uncertaintie s that need to be properly 
modeled (see e.g. lRhodes et al.lE^007k . 

Following an approach that turned o ut to be extremely suc- 
cessful in cosmic sh ear studies (see e.g 'Heyma ns et al.l 20061: 
Massey et al.l 12007 al) . we aim at overcoming these difficulties 



and facilitating the comparison between theory and observa- 
tions by performing realistic simulations, including the rele- 
vant observational eff'ects. Using suitable ray-tracing codes, the 
gravitational tensing produced by realistic mass distributions 
obtained from numerical N-body and hydrodynamical simula- 
tions can be included. Simulated images can finally be analyzed 
like real observations. 

In this paper, we describe a simulation pipeline that we will 
use to simulate observations of tensing events in galaxy clus- 
ters, but which can easily be used also for simulating tensing by 
large scale structures in wide fields. These simulations will be 
used in forthcoming papers to evaluate the accuracy of several 
tensing techniques mentioned above. 

The plan of the paper is as follows. In Sect. [2| we present 
the architecture of the simulator and describe the methods used 
to include tensing eff'ects and observational noise. In Sect. [3l 
we use a galaxy cluster obtained from N-body simulations and 
use our code to create mock optical observations of its central 



regions. We discuss how several strong tensing features can be 
extracted from the simulated images and how the properties of 
gravitational arcs can be measured. We also quantify how these 
properties depend on the assumed source models. Finally, in 
Sect m we report our conclusions. 

2. The simulator 

In this Section we describe the features included in our simu- 
lator for optical observations. In the following we will assume 
that: 

- the lenses that are eventually placed in front of the sources 
are obtained from N-body and hydrodynamical simula- 
tions. However, any deflector can be used for testing the 
codes, for example mass distributions that can be generated 
using analytic models; 

- a method for calculating the deflection angle fields out of 
the simulated mass distributions has been separately imple- 
mented. 

Our simulator is coupled to an existing ray-tracing code 



that has bee n widely used in the past (see e.g. Meneghetti et at 



I200dl200l[l2()0il2()()5b and that wifl be described in Sect.lU 
We have also implemented other routines for calculating the de- 
flection angles from a distribution of particle s using a tree algo- 
rithm (similar to that describ ed inlAubert et al.ll2007k and Fast - 
Fourier- transform methods (iBartelmann & Meneghetti '2004*). 
Ray-tracing through multiple planes for mimicking the efl'ects 
of the larg e- scale- structu re of the universe has been also imple- 
mented (P ace et a02007k . The choice of the implementation is 
not important for this paper, while it will become more relevant 
in future applications of this code. In general, for a ray-tracing 
simulation through a single plane of matter, given a deflection 
angle map a, the angular position y where a light ray is emitted 
can be easily calculated from the apparent position x through 
the lens equation: 



As ^ . X 
— a(x) . 



(1) 



where all angles are referred to an arbitrarily chosen optical 
axis passing through the observer. The angular diameter dis- 
tances between the lens and the source. As, and between the 
observer and the source, D^, have been used to re-scale the de- 
flection angles to the proper lens- source configuration. 

Assuming we know the distribution of matter along the line 
of sight and to be able to characterize its deflection field, we 
proceed now to generate a population of source galaxies, to as- 
sign them a given morphology and surface brightness distribu- 
tion, to distort them according to the tensing eff'ect produced by 
the matter between the observer and the sources and, finally, to 
visualize the results of a virtual observation including all rel- 
evant sources of noise and the convolution with a particular 
instrument. 

2.1. Galaxy positions, luminosities and redshifts 

The virtual observation comprises a user-specified field-of- 
view. This defines the size of the light-cone along which source 
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Fig. 1. Spectral-energy distributions used to characterize the 
galaxies of different spectral type. The galaxies are divided into 
four principal classes: ellipticals/SO (type 1), early spirals (type 
2), late spirals (type 3) and irregulars (type 4). The wavelengths 
are given in A, while flux is arbitrary. 

galaxies are distributed. When generating galaxy positions, we 
randomly choose the projected position of the sources on the 
sky. Their orientation is also randomly selected. Thus, we do 
not consider, at the moment, clustering of galaxies and intrin- 
sic alignments that could represent a potential source of errors 
in weak lensing measurements. We will address this issue in 
future work. 

The source distances, morphological types and intrinsic 
magnitudes are generated using observed luminosity functions 
per given redshift bin. We choose to adopt the luminosity func- 
tions derived f rom th e VIMOS VLT Deep Survey (VVDS, 
Le Fevre et al.l (l2005b ). for which difl'erent versions exist for 
fou r principal sp ectral types of galaxies up to z = 1.5. 

IZucca et alT(l2QQ6) have divided the galaxies in four types. 



corresponding to E/SO (type 1), early spiral (type 2), late spiral 
(type 3) and irregular template (type 4). This work has been 
done by matching the UBVRI magnitud es with empirical set s 
of spectral-energy-distributions (SEDs) (lArnouts et al.lll999b . 
These SEDs are shown in Fig. [T] in difl'erent colors and have 
been used to characterize the flux of the synthetic galaxies in 
our simulations. 

The luminosity function for each spectral type and redshift 
bin is of the form 



(2) 



m)dM = O.40*lnlo(lO«-4(^*-^>)""'. 
.exp(-10«-4(^*-^>yM. 

It gives th e number o f galaxies per unit magnitude and per cu- 
bic Mpc (ISchechteil ll976). The characteristic absolute mag- 
nitude M*, the characteristic density 0* and the slope a are 



determined by fitting the observed magnitude distributions of 
galaxies with the same spectral type and with spectroscopic 
redshift in the interval betw een Zmm a nd Zmax- We use the val- 
ues reported in Table 3 of Zuc ca et akl ( 2006) for the rest frame 
B-band (Johnson B filter). The galaxies used in the simula- 
tions are generated such as to reproduce the observed lumi- 
nosity functions in this reference band. The magnitudes in the 
other bands are easily calculated by using the previous SEDs, 
suitably normalized such as to assign the correct luminosity in 
the B-band, and by convolving them with the appropriate filter 
curves. 

The luminosity functions for diff'erent morphological types 
exist for galaxies up to z = 1.5. For hi gher redshift s, we adopt 
th e luminosity functio ns published by l llbert e t al.l (2005) and 
by IPaltani etaP (l2QQ7h and we extrapolate the fractions of el- 
lipticals, spirals and irregulars from those observed at lower 
redshift. In particular at z ^ 2 we assume that only 5% of the 
galaxies are ellipticals or SO, while the vast majority of them 
are spirals and irregulars. 

Certainly, important assumptions are made in the process 
of generating the population of galaxies filling the light cones. 
First, the VVDS is limited to /-magnitude 24 (AB magnitude). 
Thus, the faint tail of the distribution of galaxy luminosities is 
extrapolated from a fit performed on brighter galaxies. Second, 
the absolute magnitudes in bands diff'erent from the reference 
B-band are calculated by adopting a limited number of SEDs. 
Third, as explained earlier, we can use luminosity functions per 
difl'erent spectral types only at z < 1.5. For larger redshifts, 
we must use a single luminosity function and extrapolate from 
lower redshift the fractions of galaxies with a given SED. 

Although it is natural to expect a dominant fraction of ir- 
regulars and spirals at high redshift, our procedure may lead 
to significant errors especially when simulating deep observa- 
tions. In order to check and validate our methods, we make 
the following test. We generate a population of galaxies as de- 
scribed above and compare the number counts per magnitude 
bin and squared arcmin to those measured in the deep obser- 
vations like the HUDF. We find a good agreement between the 
number counts independent of the band selected for the experi- 
ment. For example, we show in the left panel of Fig.[2lthe com- 
parison between our prediction based on the VVDS luminosity 
functions (solid line) and the HUDF (dotted line) in the z-band 
(F850LP filter on board the Hubble telescope). Although the 
selected band is much redder than the reference band and the 
size of the HUDF is rather small, the curves are almost per- 
fectly superposed, confirming the reliability of our results de- 
spite the assumptions made. 



2.2. Galaxy surface-brightness distributions 

The light emission from galaxies falling into the telescope field 
of view is modell ed using shapelet fu nctions. This approach 



or view IS moaeii ea usmg snapeiet ru nctions. ims approacn 
has been used by iMassey et al.l ( 2QQ4I) , although us ing a dif- 
ferent shapelet decomposition method. According to lRefregieJ 
(l2QQ3h . the galaxy two-dimensional surface-brightness distri- 
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Fig. 2. Left panel: Number counts per apparent magnitude bin and squared arcmin predicted by our method (solid line) and 
observed in the HUDF (dotted line). Magnitudes are given in the z-band. Right panel: probability distributions of the galaxy 
half-light radii in the HUDF (dashed histogram) and in a field simulated with Sky lens (solid histogram). 



bution, /(jc), can be approximated by a summation over a set of 
basic functions, 



I(x)= Yj hBn{x-Xc\p) 



(3) 



where Xc indicates the galaxy centroid, x = (xi, X2) is a two- 
dimensional vector, and n = (^1,^2). The basic functions 
Bn(x', P) are called shapelets and have the form 

Bnix^P) =fi-'(f>n,(fi-'xi)(f>n,(j3-'x2) , (4) 

where the one-dimensional functions (pn(x) are related to the 
Gauss-Hermite polynomials, Hn(x), 

cPn(x) = Y2r7i^l\\r^l^Hn{x) exp(-xV2) . (5) 

The parameter p defines the scale of the galaxy image. The 
coefficients In properly weight the contribution of each shapelet 
function to the surface-brightness of the galaxy and are given 
by 



I(x)B„(x;l3) . 



(6) 



Basic photometric information about the galaxy is readily 
derived from the shapelet coefficients. For example, the total 
flux / is given by 



even / \l/2 / \l/2 

-i(2-«,-«2)( «1 1 / «2 



^ \ni/2 \«2/2 



(7) 



and the rms radius r/, defined as = J d^xx^I(x)/I, is given 
by 

even 



ni 



1/2 



1/2 



4. 



(8) 



The shapelet decomposition has pr oven to be a power- 



ful technique for image anaysis (see e.^. | Melchior et al 
' ^l2005h 



2007; 



Massev et al.l I2007bl: iGoldberg & BaconI 2005h . Among the 
many applications, Kelly & McKay ('2004") have shown that it 
can be applied to morphological classification. In fact, they find 
that galaxies of diff'erent morphol ogical types separat e well in 
shapelet space. On the other hand, Massey et al. ("2004") showed 
that an infinite number of synthetic galaxy images can be eas- 
ily created from a set of shapelet coefficients obtained from 
decomposing real galaxies. Following these ideas, we aim at 
using real shapelet decompositions to characterize galaxy mor- 
phologies in our simulator. In detail, we proceed as follows: 



we use the newly developed code by lMelchior et al.l (l2007h 
for shapelet-decomposing the ima ges of ~ 3000 galaxie s 



1. 



from the GOODS HST/ACS data (iGiavalisco et al.ll2004D . 
For each galaxy, the code provides the best fit set of co- 
efficients, /^^^^^, together with their errors, and the scale 
parameter, yg^ooDS . 

2. for about half of the galaxies in the database a spectral clas- 
sification is available, which allows us to distinguish be- 
tween morphological types and to assign a SED to each 
galaxy. For those galaxies which have no classification, we 
assume that they are irregulars or late spirals, since they are 
typically faint and distant blue objects. Based on this clas- 
sification, the galaxy decompositions are divided into four 
catalogs corresponding to the SEDs described in Sect. I2.lt 

3. aiming at generating a galaxy of type M, one entry from the 
corresponding catalog is randomly chosen. Since our cata- 
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logs contain a limited number of galaxies, which may be 
smaller than the number of galaxies to be generated in the 
simulation, replications of the same galaxy may occur. To 
avoid that, whenever a galaxy is chosen, we slightly change 
the values of the shapelet coefficients (within their errors) 
to obtain a new, slightly different, galaxy image; 

4. the galaxy image is finally rescaled on the basis of the input 
flux. From the decompositions contained in our shapelet 
database we build an empirical relation between the scale 
parameter and the flux. Then, we use this relation and its 
dispersion to draw the scale parameter (3 that corresponds 
to the input flux /. Eq. [7] however states that / oc pin. Thus, 
to ensure consistency, when scaling a GOODS galaxy with 
scale parameter /^o and flux /q to one with scale parameter 
and observed flux /, we also adjust its shapelet coefficients 
by multiplying them by (I/Io)(j3o/l3). In order to verify the 
reliability of our method, we compared the size distribu- 
tion of simulated galaxies with that observed in the HUDF. 
The probability distribution fu nctions of the half-light r adii 
(as outputted by SExtractor, iBertin & ArnoutsI ( 19961) ) in 
the simulated and in the observed HUDF are shown in the 
right panel of Fig O Fo r the observed HUDF, we refer to 
the catalog published by lBeckwith et al.l (l2006|) . The simu- 
lated images have been analyzed consistently, i.e. by adapt- 
ing the SExtractor parameters to those used in the obser- 
vations. Clearly, the distributions are very similar, although 
the half-light radii tend to be slightly larger in the simu- 
lation (median Rh = 0.211") than in the real observation 
(median Rh = 0.204"). An even more accurate result will 
be achieved when our shapelet database is extended using 
deeper observations, including a larger number of low-flux 
sources; 

5. once a set of coefficients and a scale parameter have been 
defined, they can be used to generate a new galaxy surface- 
brightness distribution. In order to ensure random orienta- 
tions, we rotate by a random angle i// e [0, tt] 

In Fig. Owe show an example illustrating the generation of syn- 
tetic galaxies out of a real GOODS galaxy. In the left panel, we 
show a spiral galaxy from GOODS. The galaxy is decomposed 
into shapelets with maximum order m = n2 = 9. The decom- 
position is then used to reconstruct the original galaxy image 
in the central panel. Finally, a second galaxy surface brightness 
distribution is generated from the same decomposition but after 
having slightly changed the values of the shapelet coefficients. 
As explained above, the coefficients and the scale parameter 
may also be rescaled to change the output flux of the sources, 
allowing us to create synthetic galaxies in arbitrary numbers. 
Moreover, as it will be shown later, using shapelets permits us 
to easily incorporate lensing efl'ects in the simulations. All im- 
ages already contain several sources of noise that will be dis- 
cussed in the following sections. 



2.3. Observing deep fields 

Once a number of galaxies have been generated, having a lumi- 
nosity, a SED, a redshift and a shape assigned, we can simulate 
observations for a particular instrumental set-up. This implies 



calculating the number of photons coming from a patch of the 
sky which are collected by each element of the CCD camera, 
after the light has passed through a sequence of mirrors and 
filters, and, eventually, through the Earth's atmosphere. 

2.3.1. Photon counts 

For the implementa tion of our simulator , we follow the pre- 
scriptions given by iGrazian et al. I (12004') in the construction 
of the Large -Binocular-Telescope Camera Image Simulator 
(LBCSIM). Given a telescope of diameter D, the number of 
photons collected by the CCD pixel at x in the exposure 
time ^exp, from a source whose surface brightness is I(x) (erg 
s"^cm"^Hz"^arcsec"^), is 



ny(x) 



4h 



J^d., 



(9) 



where p is the pixel size in arcsec, h is the Planck constant 
and T(A) is the total transmission, which is given by the prod- 
uct of the atmospheric extinction for a given airmass. A' (A) = 
^Q-o.4 airmass A(^) (airmass= for observations from space), and 
the efficiencies of the CCD, C(A), of the filter, F(A), of the mir- 
ror, M(/l), and of the optics 0(A), 



T(A) = A\A) • C(A) • F(A) • M(A) • 0(A) . 



(10) 



The number of photons is converted into a number of ADUs 
(Astronomical Data Units) by dividing by the gain g of the 
CCD: 



ADU(x) = 



ny(x) 



(11) 



2.3.2. Lensing effects 

Lensing by intervening matter along the line of sight is applied 
by mapping the CCD pixel coordinates using the lens equation 
[U The pixel at x is now connected to the position y on the plane 
where the galaxy surface brightness is assigned. Then, the sur- 
face brightess at j is converted into ADUs using Eqs.[9landfTn 
and the computed value is assigned to the pixel at position x on 
the CCD: 



ADUm(x) = ADUly(x)] , 



(12) 



Note that in the strong lensing regime one coordinate y can 
correspond to more than one pixel coordinate x. 

The application of Eq. [121 in the case of arbitrary pixel 
scales involves two steps. First, since the deflection angles may 
be determined on a grid whose point positions generally difl'er 
from those of CCD pixels, the deflection angle map must be in- 
terpolated at each pixel position, x. We use a bi-cubic interpola- 
tion scheme, which works well provided the spatial resolution 
of the virtual CCD is not too large compared to that of the input 
deflection angle map. Second, since the galaxies are distributed 
in redshift, the deflection angles must be rescaled by the factor 
As/^s appearing in Eq. [T]for each source. Especially when 
simulating observations of large fields, this operation may be- 
come time consuming because of the large number of pixels 
and sources. To shorten the calculations, we bin the sources in 
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Fig. 3. Example of generation of synthetic galaxies with shapelets. In the left panel the original image of a GOODS galaxy is 
shown. The galaxy is decomposed in 9 x 9 shapelets and reconstructed in the central panel. A second galaxy is then generated 
from the same shapelet decomposition by slightly modifying the coefficients and displayed in the right panel. The same color 
scale is used for all three images. To reproduce the background and the noise properties of the image in the left panel, we have 
simulated an observation with HST/ACS (filter F850LP, ^exp = 10600 sec). Each frame is 2.9'' on a side. 



redshift and project all the galaxies in a given redshift bin on 
a single plane. For example, we assign to all galaxies at a red- 
shift between Zmin and Zmax the redshift Zs = (Zmin + Zmax)/2. 
By changing the size of the redshift bins, one can increase or 
decrease the number of source planes. All sources falling into 
the same redshift bin can then be processed at the same time, 
saving a considerable amount of computational time. 

Note that using the shapelet formalism for representing 
the source galaxies allows to apply weak lensing transforma- 
tions through relatively simpl e operators acting on the shapelet 
states. For first order lensing. iRefregied ( 2003 ) shows that the 
transformation can be written as 

AZ)f/iensed(x) - [l + k(x)K + yi(x)S t] ADU(x) , (13) 

where k and y = (71,72) are the convergence and the shear, 
respectively. The operators K and § can be expressed in terms 
of the raising and lowering operators for the quantum harmonic 
oscillator as 



K = [af -h al^ -a\- al\ 

§ I = - ^a[^ - c^^ - al -\- ^2 j 
02 = ci^ (2^ — a^a 



(14) 

(15) 
(16) 

Thus, by using these formulas in combination with Eqs. O [9] 
and [m the lensing distortion can be expressed in terms of 
transformations of the shapelet coefficients, /„. 

2.3.3. Convolution with the PSF and seeing 

Several sources of noise can then be added to the image to re- 
semble real observational conditions. 

First, the images can be convolved with the instrumental 
point- spread-function (PSF): 



The convolution can be done in two ways. In analogy to the 
galaxy images, a given PSF model, F{x), can be decomposed 
into shapelets. Convolution is entirely analytically possible 
in shapelet space, given the invariance of the shapelet func- 
tions unde r Fourier t r ansfor mations. The full formalism is ex- 
plained in Refreded (|2003|). This approach has been used by 
Massey et al. ( 2007a ^ for image simulations. The method is ef- 



I\x) = f d^x'I(x')F(x - x') . 



(17) 



ficient for doing convolutions involving un-lensed or weakly 
lensed galaxies, i.e. galaxies with known shapelet decompo- 
sitions. For very distorted galaxies, like in the strong lensing 
regime, applying this method would require to further shapelet- 
decompose the lensed images, whose shapelet coefficients can- 
not be easily derived from the un-lensed ones (for example us- 
ing EqfT3l). Thus, this would lead to a loss of computational 
efficiency. For this reason, when simulating fields including 
strong lensing features, we opt for applying standard EFT tech- 
niques for doing convolution operations. 

This approach allows to simulate the shape of the PSF very 
realistically. However, the PSF is unique for the whole image, 
i.e. local variations of the PSF shape were not mimicked so far. 

In weak lensing applications, the PSF needs to be corrected 
in order to extract the lensing signal from the source elliptici- 
ties. Residuals in the PSF corrections, due for example to the 
above mentioned variations o f the PSF shape are one of the 
major sources of error (see e.g. 'Rhode s et al.i2007|) . 

To take this into account in our simulations, we propose the 
following approach to model the PSF. The general idea is that 
distortions of the PSF shape can be described as an additional 
"lensing" eff'ect. Anisotropics cause the images of stars to ac- 
quire irregular shapes, being stretched along some particular 
directions as if an external shear was applied. Such an external 
shear can be characterized by a power spectrum, whose ampli- 
tude defines the intensity of the anisotropics and whose shape 
characterizes the spatial scales over which it varies. 

Let the shear power spectrum be Py. We can link the power 
spectrum to a source potential power spectrum using standard 
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definitions. Since the components of the shear are combina- 
tions of the second derivatives of the lensing potential, 



1 



72 



dxidx2 

in the Fourier space we have that 



n 
n 



^1 ^ 

2 



(18) 
(19) 



(20) 
(21) 



where the hat denotes Fourier transforms and A: = (^i, ^2) is the 
wave vector. The power spectrum of the source potential is thus 



(22) 



Using this formalism, the distortions of the PSF can be in- 
troduced by applying additional lensing to the images already 
convolved with a PSF function. Given a potential if/, the corre- 
sponding deflection angle field is 



and in Fourier space 
a2 = • 



(23) 

(24) 
(25) 



Thus, the power spectra of the two components of the deflec- 
tion angle are: 



— P 

For simplicity, we assume that Py is a Gaussian, 



Pa-, — ~ - " ' 



Py(k) = Ay exp - — 



(26) 
(27) 

(28) 



whose amplitude is Ay and whose characteristic scale is ^cut- 
This last parameter sets the scales over which spatial variations 
of the PSF shape occur. 

Thus, we use Py to generate maps of the deflection angle 
components. For doing that we assume that these maps are 
Gaussian, with zero mean and variance defined by the ampli- 
tude of the input power spectrum. The deflection angles are 
used to lens the images as explained in the previous section. 

When simulating ground based observations, seeing is sim- 
ulated by further convolving the image with a Gaussian G of 
size (t: 



d^xI\x')G(x - x') . 



(29) 



In Fig. m we show some examples. In the top panel, an array 
of stars convolved with an isotropic PSF with Gaussian profile 
and FWHM of 0.6" (including seeing) is shown. The angular 



m m 






1 










m 


1 


m m 



Fig. 4. Images of an array of stars observed with an LBT- 
like telescope. The angular separation between the stars is 
10" along the x axis. In the top panel, an isotropic PSF 
(FWHM=0.6") is mimicked. In the middle panel, anisotropics 
varying on scales of ~ 1" are included. In the bottom panel we 
simulate larger- scale variations of the PSF shape on scales of 
15". 



separation of the stars along the x axis is 10". In the middle 
panel, we introduce variations of the PSF shape on scales of 
~ 1". The amplitude of the distortions is intentionally quite 
large and correspond to = 0.9. Instead, in the bottom panel 
we mimic a PSF that slowly varies on scales of ~ 15". These 
examples show that by properly choosing the parameters Ay 
and ^cut a suitable level of noise in the PSF can be added. 

2.3.4. Sky background and photon noise 

The number of ADUs per pixel from the sky background are 
given by 



4hg 



P 



-dA, 



(30) 



where S (A) is the sky flux per square arcsec at a given wave- 
length and T(A) is calculated now at airmass = 0. The sky 
flux is modeled difl'erently for observations from the ground 
and from space. For ground based observations, we use dif- 
ferent night sky spectra according to t he difl^erent phases o f 
the moon. The spectra are those used in iGrazian et al.l (2004). 
They correspond to 0, 3, 7, 10 and 14 days after the dark moon. 
For simulating the sky as seen from the space, we adopt sev- 
eral spectra of zodiacal light, representing the typical varia- 
tions of the zodiacal ba ckground radiation and of the earthshine 
(Giava Hsco et al.ll2002h . 

Photon noise is assumed to be Poissonian. For each pixel 
the noise (in ADUs) is generated by adding random deviates 
with Gaussian probability distribution of mean zero and width 



CTNix) 



^exp 



RON\^ ADU(x) -\- ADU,^y 
g j g 
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'^exp 
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(31) 
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Fig. 5. The same galaxy displayed in the central panel of Fig. [3] 
is shown here as it should be observed in the /-band by LBT 
with a diffraction limited PSF (left panel) and in bad seeing 
conditions (cr = 1.2'', right panel). The exposure time is ^exp = 
2000 sec and the airmass is 1. The scale of each frame is 2.9''. 




In the previous formula RON is the read-out noise of the chip, 
^exp is the number of exposure s, a is the flat-field te rm, which 



we fix at a = 0.005 following iGrazian et al.l (120041) . The term 
/ indicates the flat-field accuracy, which is determined by the 
number of flat-field exposures and by the level of the sky back- 
ground as 



f = (NffB-g)- 



(32) 



A more detai led explanation of these formulas is given in 
Section 4.2 of iGrazian eTaP (|2004 . 



2.3.5. Producing the images 

By properly setting the total efficiency (using the correct eflft- 
ciency curves for optics, mirrors and CCDs) observations from 
several instruments can be simulated. As an example, we show 
in Fig. [5] how the same galaxy displayed in the central panel 
of Fig. [3] would appear when observed with the LBT with a 
diff'raction limit PSF (left panel) and in bad seeing conditions 
(cr = 1.2"; right panel). The exposure time is 2000 sec with 
airmass = 1. The pixel size of the LBT camera is 0.224", in 
contrast to the 0.03" of HST/ACS (after drizzling). The spi- 
ral structure of the galaxy is thus not resolved by LBT. Seeing 
blurs the image, destroying most of the remaining morpholog- 
ical information. 

Combining several galaxies into the same field, we are able 
to produce artificial galaxy fields. One example is given in 
Fig.[6l where we simulate a deep field of 1' x 1' observed with 
the HST/ACS with an exposure time of 10600 sec. In this simu- 
lation, no lensing eff'ects have been included. As demonstrated 
here, using our procedure for generating artificial galaxies we 
can produce very realistic simulated images of the sky. 

2.3.6. Foregrounds 

The light emission from foreground galaxies can be included 
in simulated images adopting the same procedure employed 
for background sources. We note that cluster galaxies in par- 
ticular represent a potential problem for lensing analyses. For 
weak lensing studies, misidentification of cluster members as 




Fig. 6. Simulation of a deep field observed by HST/ACS in the 
z-band. The exposure time is ^exp = 10600 sec and the field of 
view is 1' X 1'. 



background sources can bias the shear measurements. In ad- 
dition the bright galaxies concentrated near the cluster centre 
(which corresponds to the strong lensing regime) complicate 
the detection of lensed sources. 

In order to include realistic galaxy clusters in our simulated 
images, we have used results from a semi-analytic model of 
galaxy formation coupled to a cluster N-body re- simulation. 
The semi-analytic niodel used in this study is described in 
IPe Lucia & Blaizo^ (l2007h and we refer to the original study 
(and references therein) for more details about the physical pro- 
cesses modelled. The model is used to generate a catalogue 
containing positions and luminosities of model galaxies at dif- 
ferent redshifts. This information is used to simulate galaxy 
images as described in the previous section. As in previous 
work, we determine the morphology of our model galaxies by 
using the B-band b ulge-to-disc ratio together with the obser- 
vational relation by ISimien & de VaucouleursI (11986) between 
this quantity and the galaxy morphological type. 

An example of a simulated observation of a galaxy cluster 
is shown in Fig. [71 The cluster is at redshift z ~ 0.3. It is ob- 
served with HST/ACS in the z-band with an exposure time of 
1000 sec. The cluster is dominated by some very bright galax- 
ies (L ~ lO^^L©) in the central part. They are expected to 
complicate significantly the observability of lensing features, 
in particular radial arcs, occurring near the cluster center. 



3. Applications of the method 

In this Section, we discuss one possible application of our sim- 
ulator, i.e. simulating gravitational arcs in the center of galaxy 
clusters. Then we apply some techniques, that can be used also 
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Fig. 7. Virtual observation in the z-band by HST/ACS of a 
galaxy cluster at z ~ 0.3. The exposure time is ^exp = 1000 
sec and the field of view is 100" x 100". The spatial distribu- 
tion and luminosities of the cluster galaxies have been obtained 
from semi-analytic models. The lensing eff'ects are disabled in 
this simulation. 

in real observations, to measure the properties of the arcs de- 
tected in the images. 



3.1. The cluster model 

The cluster used in this 



ously studied by iDolag et a 



pa per is one the clusters previ- 
" (l2005h. Its le n sing p roperties 



were investigate d also by 



Puchwein et al.l (|2005|) and by 



Meneghetti et al.l (l2007|) . A detailed description of the sim- 
ulation can be found in these papers. Its reference name is 
gl. Although several versions of this cluster exist that include 
the treatment of the gas component with several physical pro- 
cesses, we have chosen to use for the present paper the sim- 
plest simulation, where the cluster contains only dark matter 
particles. 

The halo has a mass of ~ 1.4 x 10^^ H'^Mq (resolved by 
more than a million particles within the virial region), thus it 
represents a very massive and efficient strong lens at z = 0.3. 
We have chosen this redshift because it is close to where the 
strong lensing efficiency of clusters is the largest for sources at 
Zs > 1 (Li et al.ll2005i,) . 

Although we are investing the lensing property of this halo 
only at at z = 0.3, during the simulation 92 time slices were 
saved from redshift 60 to 0. These are equispaced in time. For 
each snapshot, we computed group catalogues and their em- 
bedded substructures using a standard friends-of-friends algo- 
rithm with a linking length of 0.2 in units of the mean particle 
separation, and a mod ified version of the algorithm SUB FIND 
(|Springel et al.l l200l|) that was extended for future applica- 



Fig. 8. Surface density map of cluster g\ corresponding to pro- 
jecting the cluster along the z-axis of the simulation box. The 
side-length of the image is 1.5 /z"^Mpc comoving (correspond- 
ing to ~ 372"). 



tions to hydro- simulations. Substructure catalogues were then 
used to construct merger history trees as described in ISpringell 
(l2005h and lDe L ucia & BlaizotI (l2007|) . These merger trees are 
the basic input needed for the semi-analytic model described in 
Sect. [2321 

The surface density map corresponding to one projection 
of the cluster is shown in Fig. [H The size of the image is 
(1.5 /z"^Mpc)^. The inner region appears quite elliptical and de- 
void of large substructures. Secondary massive clumps of mat- 
ter are located at distances > 400 /z"^kpc from the cluster cen- 
ter. 



3.2. Calculations of deflection angles 

Ray-tracing simulations are carried out using the technique de- 
scribed in detail in several ear lier papers fe.g. lBartelmann etaP 
1998: M eneghetti et al.ll2000h . 

We select a cube of 3 h~^Mpc comoving side length, cen- 
tred on the halo centre and containing the high-density region 
of the cluster. The particles in this cube are used for produc- 
ing a three-dimensional density field, by interpolating their po- 
sition on a gri d of 1024^ ce lls using the Triangular Shaped 
Cloud method (IHockney & Eastwood, 1988,) . Then, we project 
the three-dimensional density field along the coordinate axes, 
obtaining three surface density maps S/j, used as lens planes in 
the following lensing simulations. 

The lensing simulations are performed by tracing a bun- 
dle of 2048 X 2048 light rays through a regular grid, covering 
the central sixteenth of the lens plane. This choice is driven by 
the necessity to study in detail the central region of the cluster. 
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where critical curves form, taking into account the contribution 
from the surrounding mass distribution to the deflection angle 
of each ray. 

Deflection angles on the ray grid are computed as follows. 
We first define a grid of 256x256 "test" rays, for each of which 
the deflection angle is calculated by directly summing the con- 
tributions from all cells on the surface density map S/^, 

where A is the area of one pixel on the surface density map 
and Xh,k and Xij are the positions on the lens plane of the "test" 
ray (/z, k) and of the surface density element (/, j). We avoid 
the divergence when the distance between a light ray and the 
density grid-point is zero by shifting the "test" ray grid by half- 
cells in both directions with respect to the grid on which the 
surface density is given. We then determine the deflection angle 
of each of the 2048 x 2048 light rays by bi-cubic interpolation 
between the four nearest test rays. 



T 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1- 




3.3. Observations 

For simulating observations of our numerical cluster, we 
consider an instrument with the characteristics of the 
planned space telescope DUNE (Dark Universe Explorer, 
|http : / /www ■ dune -mission . netj . This mission was re- 
cently proposed to ESA within its "Cosmic Visions" pro- 
gramme. It is targetted at studying the dark components of the 
universe with a wide field imager, in particular through weak 
lensing. However, thanks to the good spatial resolution, the 
panoramic field of view (0.5 square degrees) and high sensi- 
tivity, this instrument will be extremely useful for many other 
studies including strong lensing, galaxy formation and evolu- 
tion, baryonic acoustic oscillations and even planet searches. 

The motivation for adopting a future instrument for our vir- 
tual observation is twofold. On one hand, we would like to il- 
lustrate that the simulator is a powerful tool for testing tech- 
niques that are commonly applied to observational data. On the 
other hand, we would like to highlight that codes like this can 
help to define the scientific tasks that can be fulfilled by future 
missions. 

DUNE is planned to have a mirror diameter of 1.2 meters, 
a spatial resolution of 0.10"/pixel in the visible and a total effi- 
ciency close to ~ 70%. During the proposal preparation, three 
broad-band filters where used, namely 8iU-\-g-\-r (ugr), a r-h/-hz 
(riz) and a / -h z + j (izy) filter. In the proposed version, DUNE 
will not provide 3 bands in the visible. On the contrary, there 
will be one filter in the visible (riz) and 3 filters in the NIR 
(Y,J,H). The tests that are discussed below have been per- 
formed using the riz filter. Just for the purpose of illustrating 
the capabilities of the simulator for mimicking observations in 
multiple bands, we used the ugr and the izy filters (see Fig.fTOl 
below). The throughputs that we assume, including the quan- 
tum efficiency of the detector and the optical design, are shown 
in Fig.H The PSF is expected to have a FWHM of 0.23" with 
ellipticity smaller than 6%. 

A virtual observation of the inner 100" x 100" region of 
cluster ^1 is shown in Fig.[TOl It is a composite ugr, riz and izy 



Fig. 9. Throughputs assumed in carrying out the simulations 
with the DUNE satellite. Three broad-band filters are used: 
ugr (dashed lines), riz (solid lines), izy (courtesy of A. Amara, 
J. Rhodes and the DUNE collaboration members). Note that 
the figure does not reflect the final choice for the filters to be 
mounted on DUNE (see text for more details). 



image, corresponding to an exposure time of 1000 sec. In these 
tests, we assume an exposure time of 1000 seconds, although it 
has been proposed that DUNE exposures will be longer (1500 
seconds). The galaxies were drawn from our the shapelet li- 
brary, as explained earlier. In order to facilitate the formation 
of strong lensing features, we intentionally placed seven "test" 
galaxies close to the caustics of the lens. They are shown in 
the left panel of Fig. [TT] where we run a simulation without 
the foreground lens and we zoom over the region where the 
caustics are located. The caustics are also shown in blue. For 
simplicity, we assume that the "test" sources are all atZs = 1.5. 
Two galaxies, labeled A and B, are placed very close to the 
cusps of the tangential caustic, thus they are expected to pro- 
duce extended cusp arcs. Two other galaxies (F,G) are placed 
along the radial caustic and their images will be radially elon- 
gated arcs. The remaining "test" sources (C,D,E) are located 
along the fold of the tangential caustic, thus they should be 
lensed to form fold arcs. The lensed images are presented in 
the central panel of FiglTTI 

Some of the arcs, and especially the radial ones, are par- 
tially or totally hidden by the light of the foreground galaxies, 
as can be seen in Fig{TOl On the basis of the cluster merger tree, 
the semi-analytic models predict the presence of a very bright 
galaxies near the cluster center. These render accurate measure- 
ments of arcs shapes behind them problem atic. In real cases , 
these galaxies should be removed (see e.g. ISand et aDl2004l) . 
Thus, we subtract the foreground galaxies from the image by 
fitting their surface brightness profiles with Sersic profiles. This 
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Fig. 11. Left panel: positions of the sources relative to the caustics of the lens. The image is 25" x 25" wide. Central panel: 
riz image of the cluster gl after subtracting the cluster galaxies using magnitude and size cuts (see text for more details). The 
size of the image is 100" x 100". The images of the sources in the left panel are indicated by the letters, followed by numbers 
corresponding to image multiplicity. The image Gl is not visible, but we mark in red the position where it should be. Right panel: 
the points belonging to the lensed images are marked in green. 




letters followed by numbers corresponding to their multiplic- 
ity. Note the radial image labeled as Gl (one of the images of 
source G) is not detectable even after subtraction of the fore- 
ground galaxies. 



Fig. 10. Composite ugr+riz+izy image of the cluster g\. 
Several strong lensing features are present in the image. The 
FOV size is 100" x 100". 



can b e done with a suitable software like GALFIT (iPeng et al 
20021) . 



We show the image after removing the foreground light in 
the central panel of Fig.[TT] Superposed to the lensed field, we 
display the cluster critical lines. Cluster members have been 
identified by selecting galaxies with magnitude ruriz < 22 and 
rms radius r/ > 0.5". Comparing to the input catalogs of clus- 
ter members, we find that foreground galaxies were efficiently 
identified by applying these selection criteria. The images of 
the "test" galaxies displayed in the left panel are indicated by 



3.4. Arc properties 

In some applications of strong cluster lensing, such as arc 
statistics, it is important to accurately measure the properties 
of gravitational arcs. Properties like length, width and curva- 
ture radius of the i mages can be used to constrain cosmological 



parameters dBartelmann et al . 1998, 2003: Wambsganss et al 
'2004';'Oguri et al.'"2003';'Li et al. 2005: Meneghetti et al.ll2005l: 



Horesh et al. 2005; Hilbert et al. 2007) but also for dete rmining 
the inner- structure of galaxy clusters ( Sand et al ] l2005h . In this 
section we illustrate a procedure to measure the shape of arcs 
that can be applied to both real and simulated, thus favoring the 
comparison between observations and simulations. 

First, we identify the pixels in the image that belong to the 
lensed images. To do that, we make local measurements of the 
sky background and set a limit above which pixels can safely 
be assigned to the arcs. For avoiding the inclusion of noise, 
especially where foreground galaxies have been subtracted, we 
select only pixels above S/N > 5. The result of this selection is 
presented in the right panel of Fig. [TTl Several images cannot 
be classified as arcs (for example the images C2, D2, etc), thus 
we exclude them from the following analysis. We focus on the 

arcs Al, Bl, CI, Dl, El and Fl. 

(Il993h (s ee also 



As suggested by iMiralda-Escudi 



BartelmannetaDi 19981: iMeneghetti et alTl2000l 120011) . we de 



fine three characteristic points on the arcs, namely 

1 . the brightest point, that is likely to be an image of the source 
center; 

2. the point at the largest distance from 1); 

3. the point at the largest distance from 2). 
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Fig. 12. Fit of the arc Al: three characteristic points have been 
identified on the image (CPl-3), through which a circular seg- 
ment (CP2-CP3) have been traced. The curvature radius (Rc) 
and the centre of curvature (CC) of the circular segment define 
the radius and the centre of the arc. 
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Fig. 13. Transversal width profile of arc Al (red solid line). 
The blue dotted line shows an interpolation of the profile with 
a Bezier curve. The horizontal lines correspond to various de- 
terminations of the arc width. 



In the case of arc Al, these three points are shown in FigfT2l 
We trace a circle segment through these characteristic points. 
Its length and curvature radius are identified to the length and 
the curvature radius of the arc, respectively. 

Measuring the arc width is less simple. Arcs are typically 
structured, as they originate from mergers of multiple images 
of the same source. Thus, even in the case of simple, elliptical 
sources, the width is not constant along the arc. Moreover, arcs 
form in crowded regions of clusters. As shown in the previ- 
ous examples, the light from cluster members can influence the 
detectability of the lensed images or of part of them, thus af- 
fecting the measurements of the arc widths and lengths. Cluster 
galaxies can also alter the shap e of arcs by mean of their lens- 
ing effects (.Meneghetti et al.ll20QQl Last but not less important, 
photon noise causes the edges of the images to be irregular. 



The transversal width profile of arc Al is shown in Fig. [13] 
(solid line). This has been recovered by making radial scans 
of the arc along straight lines passing through the center of 
curvature CC (see Fig. [12]) and intercepting the circular seg- 
ment CP2-CP3 at equi- spaced angular separations. The angle 
is given in degrees with respect to the line CC-CP2. At each 
scan, we measure the maximal distance between arc points in- 
tercepted by the straight line. Some large scale fluctuations of 
the width can be recognized, that are caused by mergers of mul- 
tiple images. To make these large scale modes more visi ble, we 
interp olate the profile with a Bezier curve (dotted line) (|Knuthl 
119861). The profile has several intervals with positive curvature 
separated by valleys. The valleys are located at the positions 
where the arc crosses the cluster critical line, as shown in the 
central panel of Fig[TT] They roughly correspond to the regions 
where multiple images of the source have merged together. 
Additionally, many other small scale fluctuations are produced 
by noise. Thus, it appears difficult to define at which position 
the arc width has to be measured. We propose here to use the 
median width. 

In the case of observations through the Earth atmosphere, 
this effect is enhanced as both the width and the length of arcs 
are affected by the instrumental PSF and by the seeing. Since 
by definition arcs have large length-to-width ratios, the blurring 
of the images produces a much larger relative change of the 
widths than of the lengths. Thus, a correction must be applied 
to the width measurements in particular. This correction can be 
made by means of a de-convolution. Let assume that the radial 
profile of the arc is well fitted by a Gaussian, 

G(x) = Q exp - - 



2(t2 



(34) 



where x is the radial distance from the arc center of curvature. 
We define "Gaussian equivalent width" Wg given by 



(35) 



In this equation, D represents the detection limit used to select 
the arc points. 

We assume that the total PSF (instrumental PSF-h seeing) 
can be described by another Gaussian function, whose FWHM 
is equivalent to the PSF size: 

1 

GpsF = — ;= exp-— — , (36) 

ylncTpsF 2^ PSF 

FWHM = 2 ^/2ln2c^psF ^ 2354^crpsF • (37) 

We further assume that the intrinsic radial profile of the arc is 
also Gaussian: 

(x - x)^ 



Go(x) = Qo exp - 



2ctI 



(38) 



Then, the observed profile is a simple convolution of two 
Gaussians and the parameters of the intrinsic profile can be de- 
rived from the observed profile via the equations 




(To = 



cr 

Go = e— . 



(FWHM\^ 



\ 2.335 / 



(39) 
(40) 
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Fig. 14. Median radial profile of arc Al (square data points). 
Several attempts of fitting it are shown: single Gaussian (solid 
line), two Gaussians (long-dotted line) and single Gaussian fit- 
ting only the data points below half of the maximum of the 
profile (short-dotted line). 



Using Eq. [35l the de-convolved Gaussian-equivalent width 



IS 



Go 



\VG,o = 2V21n^(ro. 



(41) 



Unfortunately, the radial profile of arcs is not well fitted by 
single Gaussians. In Fig. [141 we show the median radial profile 
of the arc Al. The data points are given by squares. The solid 
line shows the best single Gaussian fit to the data. Clearly, a 
single Gaussian provides a bad fit, both in the centre and in 
the wings of the profile. In particular, in the wings the fit is 
lower than the data. This implies that the Gaussian-equivalent 
width underestimates the width of the arc. To be consistent with 
our definition of the length, this is the maximal radial exten- 
sion of the arc above a given brightness threshold. In this case, 
the threshold corresponds to ~ 100 ADUs. At that level, the 
Gaussian-equivalent width is ~ 2 pixels (i.e. ~ 0.23") less than 
the true width. 

The fit improves if a combination of Gaussians are used: 



M 



<^W = 2j2/exp- 



(42) 



The coarsely-dotted line in Fig. [14] shows the best fit obtained 
by combining M = 2 Gaussians. The reduced is smaller by a 
factor of ~ 2 in this case. The resulting fit consists of a combi- 
nation of a broad Gaussian describing the wings of the profile, 
and of a narrow one describing its central part and contributing 
very little to the wings. This suggests that, in order to mea- 
sure the width of the arc, we could use a single Gaussian to fit 
the external parts of the profile, without aff'ecting significantly 
the measurements. For example, the Gaussian that fits only the 
data points below the half of the maximum of the profile is 
shown by the finely-dotted line. Comparing to the previously 
described two-Gaussian fit, the diff'erences in the wings of the 
profile are minimal. 



Even using multiple Gaussians, the de-convolution of the 
profile can be done as described above. Indeed, each Gaussian 
component can be de-convolved individually and summed up 
to provide a de-convolved profile. The arc width can be then 
derived straightforwardly. 

The same fitting procedure can also be applied to the tan- 
gential profile of the arc, in order to correct for the PSF eff'ects 
on the arc length. However, as explained above, for arcs with 
large length-to-width ratios, the correction is much less impor- 
tant than for the arc width, since the tangential profile declines 
much more gently than the transversal profile. For example, 
the best fit Gaussians describing the decline of the tangential 
brightness profile on both sides of arc Al have cr > 15 pixels. 
The Gaussian best fitting the wings of the radial profile shown 
in Fig.[Hlhas cr = 3.4 pixels. Using Eqs. [35]and|4TI the PSF 
correction amounts to ~ 2% for the width and to < 0.3% for 
the length in that example. 

The width measurements made on arc Al are displayed 
in Fig [13] using horizontal lines. Results are shown both for 
the PSF convolved and de-convolved widths. The Gaussian- 
equivalent width obtained by fitting the arc radial profile with 
a single Gaussian is 0.99". This is ~ 15% smaller than the 
Gaussian-equivalent width measured by using two Gaussians 
for fitting the profile (1.14"). The latter is very similar to the 
true median width (1.15"). Finally, the de-convolved width, ob- 
tained from the two-Gaussians fit, is 1.12". 

The analysis made on arc Al has been repeated over the 
other arc-like images in Fig. [TT] We summarize the properties 
of all these arcs in Table [T] In column 2 we list the arc lengths 
(without PSF corrections), followed by several definitions of 
arc widths, listed in columns 3-6. We report both the maxi- 
mal and the median widths and the Gaussian-equivalent widths 
before and after PSF-correction. In column 7, we list the arc 
curvature radii. 

The comparison between the median and the maximal 
widths shows that large fluctuations of the arc thickness are 
possible along the azimuthal profiles. Diff'erences of order 30% 
are typical but the maximal width can be even 80% larger than 
the median width, as shown for the radial arc Fl (see table). 
This image is located in a noisy region, near to some cluster 
galaxies, thus the edges of the image are very irregular. This 
confirms that defining the arc width - and other arc properties 
that depend on this definition, like the length-to-width ratio - is 
a critical issue. 

The agreement between the Gaussian-equivalent and the 
median widths is generally good. The PSF-corrections are typ- 
ically of order few percent, as expected given the narrow PSF 
of DUNE. 

Some arcs are particularly straight. Two examples are the 
arcs Dl and Fl that have very large curvature radii. In such 
cases, we arbitrarily assign a curvature radius of 1000". 

By activating and/or deactivating diff'erent features of the 
simulator, the techniques used to analyze the images can be 
tested to quantify how reliable are the measurements. In the 
last column of Tab. [T] we show the median widths measured in a 
simulation where the PSF FWHM is set to zero. The remaining 
simulation parameters (instrument, exposure time, etc.) are not 
changed. Although the PSF correction is small, we notice that 
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0.76 
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0.69 


Dl 


3.91 


0.83 


0.62 


0.61 


0.59 
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0.65 


El 


8.04 


1.17 


0.93 


0.94 


0.86 


18.5 


0.73 


Fl 


4.2 


0.87 


0.48 


0.45 


0.44 


1000 


0.47 



Table 1. Summary of the main properties of the arcs shown 
in Fig. [TT] Column 1: arc names. Column 2: arc lengths. 
Column 3-6: different definitions of arc widths (maximal, me- 
dian, Gaussian-equivalent, PSF-corrected). Column 7: arc cur- 
vature radii. Column 8: arc widths in a simulation with null 
PSF. All quantities are given in arcsec. 



these newly determined widths are in a slightly better agree- 
ment with the PSF-corrected widths in column 6 (r.m.s. devia- 
tion of ~ 0.16") than with the previous median or Gaussian- 
equivalent widths in column 4 and 5 (r.m.s. ~ 0.23" and 
~ 0.25", respectively). This suggests that the de-convolution 
technique described above is correctly applied and allows to de- 
termine the true widths with an uncertainty of order one pixel. 

3.5. Detectability of arcs 

An important question that can be answered using our simula- 
tor is how do the properties of the arcs depend on source mod- 
eling. 

First of all, we wish to determine for what intrinsic fluxes of 
the sources are the arcs detectable in DUNE images, assuming 
an exposure time of 1000 seconds. For this purpose, we run 
several simulations, using the same morphological models of 
the sources, i.e. without changing their shapelet coefficients, 
their orientation and their spectral type, but changing only the 
input fluxes such that the apparent magnitudes in the riz band 
vary between m^^ = 23 and m^^ = 28. Then, we quantify if the 
arcs are detectable by counting the number of pixels belonging 
to the arcs that emerge above a minimal S/N ratio. As done 
before, we fix this limit to S /Nj^m = 5. 

The results of this test are shown in Fig. [15] for the arcs 
Al to Fl using diff'erent line-styles. In the top-left panel, we 
show the arc areas above S/Nj^m as a function of the intrin- 
sic apparent magnitude of the sources (in absence of lensing 
magnification). The limiting source magnitude for which arcs 
are detectable significantly diff'ers from arc to arc, ranging be- 
tween 24.5 and 27. The most prominent arcs, i.e. arcs Al and 
Bl, that are the most magnified sources, can be detected only 
for source intrinsic magnitudes below 25.5 and 26, respectively, 
while some shorter arcs, like arcs CI and El, can be detected 
even if the sources are fainter (rririz ^ 27). We shall recall that 
lensing does not change the surface brightness of the sources, 
but only the solid angle under which they are observed. Thus, 
the diff'erences between the detectability limits must reflect the 
difl'erent size, shape and brightness profile of the sources. 

The observed properties of the arcs, in particular the length 
and the width, obviously depend of the luminosity of the 



sources. In the top-right and in the bottom-left panels of Fig. [15] 
we show how the lengths and the widths scale as a function of 
the intrinsic magnitude. As expected, as the sources become 
fainter, the arcs shrink in both the tangential and in the ra- 
dial directions. Such decrements, however, occur at diff'erent 
rates. While the lengths decrease slowly, the change of the 
arc widths is more rapid. This is clearly due to the presence 
of several bright knots along the azimuthal profile of the arcs 
(see e.g. Fig. [13]). These originate from multiple images of the 
same source that merge from opposite sides of the critical lines. 
Thus, as long as the arcs emerge from the background, their 
length cannot become shorter than the maximal separation be- 
tween the knots. Only for very faint magnitudes do the arcs 
break into smaller pieces and the lengths decrease rapidly. This 
behavior implies that the length-to-width ratio is not a constant 
function of the source magnitude. As shown in the bottom-right 
panel, especially for the longest arcs, the length-to-width ra- 
tio is generally inversely proportional to the source luminosity, 
once the source morphology has been fixed. For example, be- 
tween niriz = 23 and m^^ = 25, the length-to- width ratio of the 
arc Al becomes larger by more than a factor of two. For arc 
Bl, it grows by ~ 65% between m^^ = 23 and rririz - 25.5. 
Also the length-to- width ratios of the arc CI and of the radial 
arc Fl grow significantly. 

To better explain the dependence of the arc appearance 
on the intrinsic properties of the sources, we run new simu- 
lations, keeping the sources at the same positions but exchang- 
ing among them the sets of shapelet coefficients that determine 
their morphology and surface brightness distributions. In other 
words, we re- simulate each arc using seven different source 
shapes. The intrinsic magnitude of the sources was fixed to 
rririz - 24. This value has been chosen such that all arcs Al-Fl 
are detectable. In the upper graph of Fig. [16] we show the area 
of each arc above 5" /A/^ = 5 as a function of the source model. 
Each model is identified by a number running from 1 to 7. The 
corresponding un-lensed images of the sources are shown at 
the center of the frames below the graphs. Note that these are 
cut-outs from the left panel of Fig. [TT] As anticipated earlier, 
by varying the morphology of the sources, different amounts of 
pixels belonging to the arcs surpass the detectability thresholds. 
In particular, the spatial extent of the arc is determined by the 
size of the sources, i.e. by their surface brightness. We notice 
that the curves corresponding to different arcs have many fea- 
tures in common. Indeed, they have several local minima and 
maxima that arise for the same source models. For example, 
all arcs have the largest sizes if the source model number 2 or 
number 7 are used, while their size is the smallest for the source 
model number 4. Looking at the corresponding source repre- 
sentations, we see that the former are low surface-brightness 
and spatially extended galaxies, while the latter corresponds to 
the most compact among the sources. 

We now consider the correlation between the source model 
and the length-to- width ratio of their images. This is shown in 
the bottom graph of Fig. [16] Compared to the previous plot, 
we notice the opposite correlation for most of the images: for a 
fixed total luminosity of the sources, the most extended galax- 
ies (the sources 2 and 7 in Fig. [16]) correspond to the lowest arc 
length-to- width ratios. 
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Fig. 15. Arc properties as a function of the source intrinsic magnitude in the riz band. Shown are the area above 
signal-to-noise ratio S /N = 5 (upper-left panel), the length (upper-right panel), the width (bottom-left) and the len 
ratio (bottom-right panel). The different line-styles correspond to the arcs Al-Fl. 
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We can interpret these results as follows. Large, low surface 
brightness images average magnification over a larger area. 
Thus, they are less distorted. Indeed, the relative changes of the 
lengths and of the widths of the images are caused by the con- 
volution of the tangential and of the radial magnifications with 
source surface brightness distribution. The distortion eff'ects are 
diluted if the source is broader. For this reason, arcs originat- 
ing from compact, faint or low surface brightness sources that 
are barely detectable are preferentially characterized by large 
length-to- width ratios. These results show that comparing the- 
ory to observations, we need to know the source population the 
arcs are drawn. As surface brightness is unchanged, we suggest 
basing the detectability criterion of arcs upon it. 

4. Conclusions 

We presented a new code for mimicking optical observations 
of distant galaxies. It includes all relevant sources of noise af- 
fecting real observations, such as sky background, photon noise 
and eff'ects of the atmosphere and of the instrumental PSF. The 
code also incorporates ray-tracing routines that allow us to in- 
clude lensing eff'ects by matter distributed between the observer 
and the sources. 



In this work we have detailed the thorough method that we 
have developed in order to model the source galaxies realisti- 
cally. The source morphologies are modeled via the shapelet 
decomposition of real galaxy images observed with HST. Four 
spectral templates are used to simulate the spectral energy dis- 
tributions of the galaxies. Their redshift and luminosity distri- 
butions as a function of their spectral classification are derived 
from observed galaxy luminosity functions in the VVDS. 

Being particularly suitable for simulating lensing by galaxy 
clusters, the code also simulates the emission from the cluster 
galaxy population, whose properties are obtained from semi- 
analytic models of galaxy formation and evolution. 

Many features of the simulator have been discussed in de- 
tail. In particular, we have shown that the code can be easily 
used to mimic observations with a variety of existing instru- 
ments, both from the space (HST) and from the ground (LBT). 

In order to illustrate the code and to introduce some poten- 
tial applications, we have simulated observations of a galaxy 
cluster obtained from N-body simulations. These virtual obser- 
vations were made using a telescope with the characteristics of 
the proposed space telescope DUNE. 
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Fig. 16. Top panel: Areas of the arcs Al-Fl above S/N = 5, 
assuming different sets of shapelet coefficients for modeling 
the sources. The source models (shown in the figures below the 
graphs) are numbered from 1 to 7. All sources have intrinsic 
magnitude ntriz = 24. Bottom panel: length-to- width ratio as a 
function of the source model 



The simulated images were analyzed using standard ob- 
servational techniques. After subtracting the light of the fore- 
ground galaxies, multiple images of several distant sources that 
were strongly lensed by the foreground cluster can be identi- 
fied. Each detected arc-like image has been classified, measur- 
ing its length, width and curvature radius. 

We discussed how the properties of arcs can be measured 
consistently in simulations and observations, in order to facil- 
itate the comparison with theoretical predictions. In particular, 
we focused on the determination of the arc widths and on how 
they can be corrected for the PSF broadening. For this purpose, 
we fit the radial profile of the arcs with multiple Gaussians. By 
comparing simulations with and without PSF convolution, we 
verified that the de-convolution technique allows one to deter- 
mine the true width of the arcs with a typical error of one pixel. 

By varying the characteristics of the sources, we studied 
the detectability limits of arcs in DUNE observations. With an 



exposure time of 1500 seconds, we expect that DUNE will be 
able to observe arcs arising from sources whose magnitudes 
are fainter than ~ 27 in the riz band. The shape of the arcs are 
found to be very sensitive to the properties of the sources. In 
particular, we found that arcs tend to acquire larger length-to- 
width ratios as their sources become fainter or more compact. 

These results are particularly intriguing for arc statistics. 
Indeed, they show that observational eff'ects may have a large 
impact on the abundance of arcs with large length-to width ra- 
tios, and they deserve careful investigation. We will dedicate a 
forthcoming paper to this subject (Meneghetti et al., in prep.). 
The source morphologies, and in particular sub- structures of 
star-forming regions, are clearly more relevant for arc statis- 
tics than is usually assumed. This is as a result of the small 
angular size of the bright star forming regions which are more 
easily distorted. Moreover, the number of very elongated im- 
ages is expected to grow rapidly as a function of the depth of 
the observations. Thus, space missions like DUNE, capable of 
making deep and wide surveys thanks to their large sensitiv- 
ity and field of view, are likely to discover a significant num- 
ber of new gravitational arcs. The usage of efficient software 
for automatic arc identification (see e.g. Lenzen et al. 2005; 
Seidel & BartelmannI I2OO7I: lAlardI I2OO6I: iHoresh et al.l 120051: 
Cabanac et al.l2007|) will facilitate the detection of strong lens- 
ing events in these large data sets. 

Acknowledgements. The A/^-body simulations were performed at the 
"Centro Interuniversitario del Nord-Est per il Calcolo Elettronico" 
(CINECA, Bologna), with CPU time assigned under an INAF- 
CINECA grant. This work has been supported by the Vigoni pro- 
gramme of the German Academic Exchange Service (DAAD) and 
Conference of Italian University Rectors (CRUI). We acknoweldge 
financial contribution from contracts ASI-INAF 1/023/05/0, ASI- 
INAF 1/088/06/0 and INFN PD5L CH acknowledges the support 
of a European Commission Programme 6th framework Marie Curie 
Outgoing International Fellowship, under contract MOlF-CT-2006- 
21891. MM thanks M. Bolzonella, G Rodighiero, E. Zucca, S. 
Bardelli, E. Vanzella and G Zamorani for helpful discussions. We 
thank A. Refregier, J. Rhodes and the whole DUNE collaboration 
for allowing us to test our code during the preparation of the DUNE 
proposal. We thank M. Barden and S. Koposov, the GOODS and 
COMB 0-1 7 team for providing the GOODS galaxy images used for 
the presented work. We are grateful to C. Peng for providing the soft- 
ware GALFIT and to C. Wolf and K. Meisenheimer for giving us the 
SEDs of the GOODS galaxies. We warmly thank V. Springel for pro- 
viding us with post-processing software. MM thanks Paolo for his co- 
operation during the preparation of the manuscript. 



References 

Alard, C. 2006, ArXiv Astrophysics e-prints 

Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 
540 

Aubert, D., Amara, A., & Metcalf, R. B. 2007, MNRAS, 376, 113 
Bartelmann, M., Huss, A., Colberg, J., Jenkins, A., & Pearce, F. 1998, 
A&A, 330, 1 

Bartelmann, M. & Meneghetti, M. 2004, A&A, 418, 413 
Bartelmann, M., Meneghetti, M., Perrotta, F, Baccigalupi, C, & 

Moscardini, L. 2003, A&A, 409, 449 
Bartelmann, M. & Schneider, P. 2001, Physics Reports, 340, 291 



M. Meneghetti et al.: Realistic simulations of lensing by galaxy clusters 



17 



Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 

132, 1729 

Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393 

Bradac, M., Schneider, P., Lombardi, M., & Erben, T. 2004, ArXiv 

Astrophysics e-prints 
Cabanac, R. A., Alard, C, Dantel-Fort, M., et al. 2007, A&A, 461, 

813 

Cacciato, M., Bartelmann, M., Meneghetti, M., & Moscardini, L. 

2006, A&A, 458, 349 
Clowe, D., Gonzalez, A., & Markevitch, M. 2004, ApJ, 604, 596 
Comerford, J. M., Meneghetti, M., Bartelmann, M., & Schirmer, M. 

2006, ApJ, 642, 39 
Dalai, N., Hennawi, J. R, & Bode, R 2005, ApJ, 622, 99 
De Lucia, G. & Blaizot, J. 2007, MNRAS, 375, 2 
Diego, J. M., Tegmark, M., Protopapas, R, & Sandvik, H. B. 2007, 

MNRAS, 375, 958 
Dolag, K., Vazza, R, Brunetti, G., & Tormen, G. 2005, MNRAS, 364, 

753 

Ellis, R., Santos, M. R., Kneib, J.-R, & Kuijken, K. 2001, ApJ, 560, 
L119 

Fedeli, C, Meneghetti, M., Bartelmann, M., Dolag, K., & Moscardini, 

L. 2006, A&A, 447, 419 
Fort, B., Prieur, J. L., Mathez, G., Mellier, Y., & Soucail, G. 1988, 

A&A, 200, L17 

Giavalisco, M., Ferguson, H. C, Koekemoer, A. M., et al. 2004, ApJ, 
600, L93 

Giavahsco, M., Sahu, K., & Bohhn, R. 2002, in HST Instrument 

Science Report WFC3 2002-12, 1 
Gladders, M., Hoekstra, H., Yee, H., Hall, R, & Barrientos, L. 2003, 

ApJ, 593, 48 

Goldberg, D. M. & Bacon, D. J. 2005, ApJ, 619, 741 
Golse, G., Kneib, J.-R, & Soucail, G. 2002, A&A, 387, 788 
Grazian, A., Fontana, A., De Santis, C, et al. 2004, PASP, 116, 750 
Heymans, C, Van Waerbeke, L., Bacon, D., et al. 2006, MNRAS, 368, 
1323 

Hilbert, S., White, S. D. M., Hartlap, J., & Schneider, R 2007, 
MNRAS, 927 

Hockney, R. & Eastwood, J. 1988, Computer simulation using parti- 
cles (Bristol: Hilger, 1988) 

Horesh, A., Ofek, E. O., Maoz, D., et al. 2005, ApJ, 633, 768 

Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863 

Kelly, B. C. & McKay, T. A. 2004, AJ, 127, 625 

Kneib, J.-R, Ellis, R. S., Santos, M. R., & Richard, J. 2004, ApJ, 607, 
697 

Knuth, D. 1986, Metafont: the Program (Addison- Wesley) 

Le Fevre, O., Vettolani, G., Garilli, B., et al. 2005, A&A, 439, 845 

Lenzen, R, Scherzer, O., & Schindler, S. 2005, A&A, 443, 1087 

Li, G., Zhang, P., & Chen, X. 2007, ArXiv Astrophysics e-prints 

Li, G.-L., Mao, S., Jing, Y. R, et al. 2005, ApJ, 635, 795 

Luppino, G., Gioia, L, Hammer, R, Le Fevre, O., & Annis, J. 1999, 

A&AS, 136, 117 
Lynds, R. & Petrosian, V. 1989, ApJ, 336, 1 
Massey, R., Heymans, C, Berge, J., et al. 2007a, MNRAS, 376, 13 
Massey, R., Refregier, A., Conselice, C. J., David, J., & Bacon, J. 

2004, MNRAS, 348, 214 
Massey, R., Rowe, B., Refregier, A., Bacon, D. J., & Berge, J. 2007b, 

MNRAS, 649 

Melchior, P., Meneghetti, M., & Bartelmann, M. 2007, A&A, 463, 
1215 

Meneghetti, M., Argazzi, R., Pace, R, et al. 2007, A&A, 461, 25 
Meneghetti, M., Bartelmann, M., Dolag, K., et al. 2005, A&A, 442, 
413 



Meneghetti, M., Bartelmann, M., Jenkins, A., & Frenk, C. 2005a, 

ArXiv Astrophysics e-prints 
Meneghetti, M., Bartelmann, M., & Moscardini, L, 2003, MNRAS, 

346, 67 

Meneghetti, M., Bolzonella, M., Bartelmann, M., Moscardini, L., & 

Tormen, G. 2000, MNRAS, 314, 338 
Meneghetti, M., Jain, B., Bartelmann, M., & Dolag, K. 2005b, 

MNRAS, 362, 1301 
Meneghetti, M., Yoshida, N., Bartelmann, M., et al. 2001, MNRAS,, 

325, 435 

Miralda-Escude, J. 1993, ApJ, 403, 497 

Oguri, M., Lee, J., & Suto, Y. 2003, ApJ, 599, 7 

Pace, R, Maturi, M., Meneghetti, M., et al. 2007, A&A, 471, 731 

Paltani, S., Le Fevre, O., Ilbert, O., et al. 2007, A&A, 463, 873 

Paulin-Henriksson, S., Antonuccio-Delogu, V., Haines, C. P., et al. 

2007, A&A, 467, 427 
Peng, C. Y, Ho, L. C, Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266 
Puchwein, E., Bartelmann, M., Dolag, K., & Meneghetti, M. 2005, 

A&A, 442, 405 
Refregier, A. 2003, MNRAS, 338, 35 

Rhodes, J. D., Massey, R., Albert, J., et al. 2007, ArXiv Astrophysics 
e-prints 

Sand, D. J., Treu, T., Ellis, R. S., & Smith, G. R 2005, ApJ, 627, 32 
Sand, D. J., Treu, T., Smith, G. R, & Ellis, R. S. 2004, ApJ, 604, 88 
Schaerer, D. & Pello, R. 2003, in Bulletin of the American 

Astronomical Society, Vol. 35, Bulletin of the American 

Astronomical Society, 1308-+ 
Schechter, R 1976, ApJ, 203, 297 
Seidel, G. & Bartelmann, M. 2007, A&A, 472, 341 
Simien, R & de Vaucouleurs, G. 1986, ApJ, 302, 564 
Springel, V. 2005, MNRAS, 364, 1105 

Springel, V, White, S. D. M., Tormen, G., & Kauffmann, G. 2001, 

MNRAS, 328, 726 
Starck, J.-L., Pires, S., & Refregier, A. 2006, A&A, 451, 1139 
Torri, E., Meneghetti, M., Bartelmann, M., et al. 2004, MNRAS, 349, 

476 

Wambsganss, J., Bode, R, & Ostriker, J. 2004, ApJL, 606, 93 

Zaritsky, D. & Gonzalez, A. 2003, ApJ, 584, 691 

Zucca, E., Ilbert, O., BardelH, S., et al. 2006, A&A, 455, 879 



